#include <iostream>
#include <stdio.h>
#define GL_GLEXT_PROTOTYPES

#ifdef __APPLE__
#include <OpenGL/OpenGL.h>
#else
#include <GL/glut.h>
#endif

#include <cuda.h>
#include <cuda_gl_interop.h>

#include <GL/glfw.h>

using namespace std;


GLuint buffer;
unsigned *pixels;
GLuint tex;
int width=512, height=512;

bool inited=false;

unsigned __constant__ lookup[256];


#  define CU_SAFE_CALL_NO_SYNC( call ) {                                     \
    cudaError_t err = call;                                                  \
    if( cudaSuccess != err) {                                               \
        fprintf(stderr, "Cuda driver error %x in file '%s' in line %i.\n",   \
                err, __FILE__, __LINE__ );                                   \
        exit(EXIT_FAILURE);                                                  \
    } }

#  define CU_SAFE_CALL( call )       CU_SAFE_CALL_NO_SYNC(call);


struct vec2 {
  float x, y;
};



typedef union colorTag {
    struct {
	unsigned char r;
	unsigned char g;
	unsigned char b;
	unsigned char a;
    };
    unsigned value;
} color;


void makeCheckImage()
{
	int i, j, c;
	static int ic = 0;  
	for (i = 0; i < height; i++) {
		for (j = 0; j < width; j++) {
			c = ((((i&0x8)==0)^((j&0x8))==0))*255;
			color cl;
			if (ic)
				c = 255 - c;
			cl.r = (GLubyte) c;
			cl.g = (GLubyte) c;
			cl.b = (GLubyte) c;
			cl.a = (GLubyte) 255;
			pixels[width * j + i] = cl.value;
		}
	}
	ic = 1 - ic;
}


void setPalette()
{
	color palette[256];

	for (int i = 0; i < 256; i++) {
		palette[i].r = i;
		palette[i].g = i;
		palette[i].b = i;
		palette[i].a = 255;
	}

	CU_SAFE_CALL_NO_SYNC( cudaMemcpyToSymbol(
		lookup,
		palette,
		256 * sizeof(unsigned)) );
}


void create_buffer(GLuint* buffer)
{
	pixels = (unsigned*)malloc(width*height*sizeof(unsigned));
}


void destroy_buffer(GLuint* buffer)
{
	free(pixels);
}

void create_texture(GLuint* tex)
{
	glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
    glGenTextures(1, tex);
    glBindTexture(GL_TEXTURE_2D, *tex);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
//    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, width, height, 0, GL_RGBA, GL_UNSIGNED_BYTE, NULL);
    glBindTexture(GL_TEXTURE_2D, 0);
}

void destroy_texture(GLuint* tex)
{
    glBindTexture(GL_TEXTURE_2D, *tex);
    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, 0, 0, 0, GL_RGBA, GL_UNSIGNED_BYTE, NULL);
    glBindTexture(GL_TEXTURE_2D, 0);
    glDeleteTextures(1, tex);
}


void display0()
{		
	static int blue = 0;
	
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glColor3f(1.0,1.0, blue * 1.0);
	glLoadIdentity();
	glBegin(GL_QUADS);
	glVertex2f(-1.0f, -1.0f);
	glVertex2f(1.0f, -1.0f);
	glVertex2f(1.0f, 1.0f);
	glVertex2f(-1.0f, 1.0f);
	glEnd();
	glfwSwapBuffers();	
	blue = 1 - blue;
}


void display()
{	
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);	
   	glEnable(GL_TEXTURE_2D);
   	glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_DECAL);
	glBindTexture(GL_TEXTURE_2D, tex);

	glBegin(GL_QUADS);
	glTexCoord2f(0.0f, 0.0f); glVertex2f(-1.0f, -1.0f);
	glTexCoord2f(1.0f, 0.0f); glVertex2f(1.0f, -1.0f);
	glTexCoord2f(1.0f, 1.0f); glVertex2f(1.0f, 1.0f);
	glTexCoord2f(0.0f, 1.0f); glVertex2f(-1.0f, 1.0f);
	glEnd();	
	glfwSwapBuffers();
}



void initGLFW()
{
	if( !glfwInit() ) {
		exit( EXIT_FAILURE );
	}
	if (!glfwOpenWindow( width, height, 0,0,0,0,0,0, GLFW_WINDOW ) ) {
		glfwTerminate(); 
		exit( EXIT_FAILURE );
	}
	
	glfwSetWindowTitle("Concentracion GLFW");
	
	glViewport(0, 0, width, height);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	glOrtho(-1, 1, -1, 1, -1, 1);
	glMatrixMode(GL_MODELVIEW);
    glEnable(GL_TEXTURE_2D);
    setPalette();
    create_buffer(&buffer);
    create_texture(&tex);
//    CU_SAFE_CALL(cudaGLRegisterBufferObject(buffer));
		
	inited=true;
}


void closeGLFW()
{
//	CU_SAFE_CALL(cudaGLUnregisterBufferObject(buffer));
	destroy_texture(&tex);
	destroy_buffer(&buffer);
	glfwTerminate();
}


__global__ void concentracion(float *c, vec2 *linec, float *s, float *tao, int line_size, float deltaL0, int Nxf)
{
  int bx = blockIdx.x;	
  int by = blockIdx.y;
  int tx = threadIdx.x;	
  int ty = threadIdx.y;	
  
  int i = bx*blockDim.x + tx;
  int j = by*blockDim.y + ty;
  
  int Nyf = Nxf;

  float dxf = 1.0 / (Nxf-1);
  float dyf = 1.0 / (Nyf-1);

  vec2 x_im1_n = linec[0]; 

  c[Nxf * j + i] = 0;

/* */
  for (int k=1; k < line_size; k++) {	    	  
    vec2 x_i_n = linec[k];
    vec2 sigma, coord;
    
    sigma.x = x_i_n.x - x_im1_n.x; 
    sigma.y = x_i_n.y - x_im1_n.y;
    float normSig = sqrt ( sigma.x * sigma.x + sigma.y * sigma.y);
    sigma.x /= normSig;
    sigma.y /= normSig;
    coord.x = dxf * i;
    coord.y = dyf * j;
    float midx = 0.5 * ( x_i_n.x + x_im1_n.x ); 
    float midy = 0.5 * ( x_i_n.y + x_im1_n.y );	  
    float aux1 = (coord.x - midx) * sigma.x + (coord.y - midy) * sigma.y;
    float aux2 = (coord.x - midx) * (-sigma.y) + (coord.y - midy) * sigma.x;
    aux1 = aux1 * aux1 / (deltaL0 * deltaL0);
    aux2 = aux2 * aux2 / ( s[k]*s[k] * (1.0 + 4.0 * tao[k]) );
    c[Nxf * j + i] += exp(-aux1-aux2) / (1.7726 * sqrt( 1.0 + 4.0 * tao[k]));
    x_im1_n = linec[k];
  }
//	unsigned result = (unsigned)(255.0 * (c[Nxf * j + i]- 6.0)/6.0);
//  p[Nxf * j + i] = lookup[250];

/*  */
}	

// Variables _h ya están inicializadas en host
int calcula_concentracion(float *c_h, vec2 *linec_h, float *s_h, float *tao_h, int line_size, float deltaL0, int grid)
{
	int size = line_size*sizeof(float);
	vec2 *linec;
	float *s, *tao, *c;

	int Nxf = 16*grid;
	int Nyf = Nxf;	
	int sizeC = Nxf * Nyf * sizeof(float);

	cudaMalloc((void**) &linec, 2*size);
	cudaMalloc((void**) &s, size);
	cudaMalloc((void**) &tao, size);
	cudaMalloc((void**) &c, sizeC);
	cudaMemcpy(linec, linec_h, 2*size, cudaMemcpyHostToDevice);
	cudaMemcpy(s, s_h, size, cudaMemcpyHostToDevice);
	cudaMemcpy(tao, tao_h, size, cudaMemcpyHostToDevice);
	dim3 DimGrid(grid, grid, 1);
	dim3 DimBlock(16, 16, 1);		
	if (!inited) {
		initGLFW();
	//	display0();
	}
	
	concentracion<<<DimGrid, DimBlock>>>(c, linec, s, tao, line_size, deltaL0, Nxf);
	cudaMemcpy(c_h, c, sizeC, cudaMemcpyDeviceToHost ); 
	float min = 1e10, max = -min;
	for (int j = 0; j < Nyf; j++) {
		for(int i  = 0; i < Nxf; i++) {
	   		float f = c_h[i*Nxf+j];
			if (f < min)
				min = f;
			if (f > max)
				max = f;
		}
	}
	printf("maxmin %g %g\n", min, max);
	for (int j = 0; j < Nyf; j++) {
		for(int i  = 0; i < Nxf; i++) {
	   		float f = c_h[i*Nxf+j];
			unsigned k = (unsigned)(255*(f - min)/(max - min));
			color c;
			c.r = c.g = c.b = 255 - k;
			c.a = 255;
			pixels[Nxf * i + j] = c.value;
		}
	}
//	makeCheckImage();
	glBindTexture(GL_TEXTURE_2D, tex);
	glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, width, 
		height, 0, GL_RGBA, GL_UNSIGNED_BYTE, 
		pixels);
	display();

	cudaFree(linec);
	cudaFree(s);
	cudaFree(tao);
	cudaFree(c);

	return 0;
}

